SEARCH TECHNIQUES FOR ROOT-UNITARY POLYNOMIALS 



KIRAN S. KEDLAYA 

Abstract. We give an anecdotal discussion of the problem of searching for polynomials 
with all roots on the unit circle, whose coefficients are rational numbers subject to certain 
congruence conditions. We illustrate with an example from a calculation in p-adic cohomol- 
ogy made by Abbott, Kedlaya, and Roe, in which we recover the zeta function of a surface 
over a finite field. 



Introduction 

In this note, we give an anecdotal discussion of the problem of searching for polynomials 
with roots on a prescribed circle whose coefficients are rational numbers subject to certain 
congruence conditions. We were led to this problem by the use of p-adic cohomology to 
compute zeta functions of varieties over finite fields; in that context, one is looking for 
certain Weil polynomials (monic integer polynomials with complex roots all on a circle of 
radius p 1 ^ 2 , for some prime number p and some nonnegative integer £), and the cohomology 
calculation imposes congruence conditions on the coefficients. In fact, the main purpose of 
this note is to show that in a particular example from [1] , the conditions obtained from the 
cohomology calculation indeed suffice to uniquely determine the zeta function being sought. 
We also illustrate with a larger example provided by Alan Lauder. 

1. Definitions 

A polynomial P(z) = YH=$ a i zl ^ ^-\ z \ °f degree n is self-inversive if there exists mgC 
with \u\ = 1 such that 

(1.1) a,i = ua~ (i = 0, . . . ,n); 

it is equivalent to require the roots of P to be invariant, as a multiset, under inversion 
through the unit circle. It appears that the class of self-inversive polynomials first occurs 
in a theorem of Cohn [3J, [5] ; it occurs naturally in the study of the locations of roots of 
polynomials and their derivatives, as in the Schur-Cohn-Marden method [131 P- 150]. 

We will call a polynomial root-unitary if its roots all lie on the unit circle; this clearly 
implies self-inversivity. This class of polynomials has been widely studied, but does not 
seem to have a standard name: the term "unimodular polynomial" refers to a polynomial 
whose coefficients lie on the unit circle, while "unitary polynomial" is often read as a synonym 
for "monic polynomial" (particularly by speakers of French, in which a monic polynomial is 
standardly a "polynome unitaire"). 
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Let P(z) e R[z] be a real root-unitary polynomial; then (11.11) must hold with either u = +1 
or u = —1, in which case we say P is reciprocal or antireciprocal, respectively. (The terms 
palindromic and antipalindromic are also sometimes used.) If deg(P) is odd, then P must be 
divisible by z + 1 or z — 1, depending on whether P is reciprocal or antireciprocal. If deg(P) 
is even and P is antireciprocal, then P must be divisible by (z + l)(z — 1). This allows 
reduction of many questions about real root-unitary polynomials to the reciprocal case. 



2. The basic problem 

The basic problem is to identify rational polynomials with roots on a prescribed circle, 
given a few initial coefficients and a congruence condition on the remaining coefficients. One 
can renormalize in order to talk about root-unitary polynomials; as noted above, there is no 
real harm in only looking at reciprocal root-unitary polynomials. In any case, here is the 
precisely formulated question we will consider. 

Problem 2.1. Fix positive integers n, k, q with n > k. Also fix positive integers mo, . . . , vn2 n 
such that rrij divides rrii for < % < j < n, and rri{ = rri2n-i for < i < 2n. Given integers 
a , . . . , a 2 „ with = q n ~ l a 2 n-i for < i < 2n, and a 2n ^ ; find all polynomials P(z) with 
all roots on the circle \z\ = y/q of the form 

2n 

P{z) = y~](ai + Cimi)z\ 

i=0 

where the Ci G Z must satisfy q = q n ~ l C2 n -i for < % < 2n, and Q = for i > 2n — k. 
The fact that Problem 12.11 is a finite problem follows easily from the estimates 

\ai + Cimi\< {^ r ^jq %l2 \a 2n \ (i = 0, . . . , 2n); 

when n is small, these estimates carry most of the information from the condition that 
P(zy/q) must be root-unitary. However, for n large, this is quite far from true. Indeed, by 
pO, Proposition 2.2.1] (see also [16J for a generalization in the context of Mahler measures), 
the space of monic root-unitary reciprocal polynomials of degree 2n has volume 

9 n n ( 2i \ n+1 ~i 9™ n o(n 2 +3n)/2 
1 I I J I ^ 1 I on+l-j 



n! - LJ - V 27 — 1 / n! - LJ - n! 

3 =l V ■> ' 3=1 

whereas the space of monic reciprocal polynomials of degree 2n whose coefficient of z l has 
norm < ( 2 ™) for i = 0, . . . , 2n — 1 has volume 

n<; ! )=n2n^sn2n^2«»"«»" 2 . 

3=1 V J ' 3=1 i=0 J j=l i=0 

For n large, these are wildly discrepant, so one expects the restriction of root-unitarity to 
carry much more information than the simple bound on the size of coefficients. 
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3. Exhaustion over a tree 



We now describe our basic approach to Problem 12.11 starting with a change of variable 
also used in [TJ. Define a polynomial Q(z) G 1i[z] of degree n by the formula 

P(z) = z n Q{z + q/z). 

Then for i — 0, . . . , n, the coefficients of z n ~\ . . . , z n of Q are obtained from a,2 n -i, ■ ■ ■ , &2n by 
an invertible linear transformation over Z. We can thus reformulate Problem 12.11 as follows. 

Problem 3.1. Fix positive integers n, k, q with n > k. Also fix positive integers mo, . . . , m n 
such that rrij divides for < i < j < n. Given integers b , . . . ,b n with b n ^ 0, find all 
polynomials Q(z) G Z[#] with all roots real and lying in the interval [— 2^/q, 2y/q], such that 

n 

Q{z) = + d t mi)z l 

8=0 

for some di G Z with di = for i > n — k. 

Our approach to Problem 13.11 is via enumeration of a certain rooted tree. 

Proposition 3.2. Fix notation as in Problem IJ.il Then there exist sets Sj C Z J for j = 

0. . . . , n — k satisfying the following conditions. 

(a) The set Sq consists of the empty {) -tuple. 

(b) For < j < n - k, if (d n -k-i, d n - k -j) e Sj, then (d n - k -i, d n - k - j+1 ) G Sj-i. 

(c) For < j < n-k, «/(d ri _ fc _i, . . . , rf n _ A; _ i ) G 5 i; then Q (z) = J27=o b i zl +Y.™=n~Lj d^z 1 
has the property that k has all roots in [—2y/q,2y/q\. 

(d) Every tuple (d n -k-i, ■ ■ ■ ,do) ^ Z n_fc such that J27=o + Y^Zo' 1 dirriiZ 1 has all roots 
in [—2^/q,2y/q] belongs to S n - k . 

Proof. Create S n - k by taking all solutions of Problem 13.11 then let Sj be the set of initial 
segments of length j occurring among elements of S n _ k . Property (c) holds by Rolle's 
theorem. □ 

We may identify a system of sets as in Proposition 13.21 with a rooted tree, where the 
children of a j-tuple in Sj are its extensions to a (j + l)-tuple in Sj + %. To solve Problem 13.11 
in practice, we perform a depth-first enumeration of such a tree, and read off the solutions of 
Problem 13 .11 as the elements of S n _ k . To describe such a tree and its enumeration, it suffices 
to describe how to compute the list of children of a given node. (One could also perform 
a breadth-first exhaustion, but in practice this seems to be inferior because of increased 
overhead.) 

Note that if one wishes to decide as soon as possible whether the number of solutions is 0, 

1, or more than 1, it may be advantageous to visit the children of a given node in "inside-out 
order" rather than in ascending or descending order. For instance, if a given tuple can be 
extended by 5, 6, 7, 8, 9, we would visit these extensions in the order 7, 6, 8, 5, 9. 

4. First approach: root-finding 

We now describe our first algorithmic approach to Problem 13.11 and its implementation 
[T2] in the case where q — 1. (See Section [8] for comments on the remaining cases.) The 
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implementation, based on an undergraduate research project with Andre Wibisono, uses the 
computer algebra system SAGE [15], and relies in particular on the following components. 

• We use the PARI C library [T3] for polynomial manipulations over Q, including 
the function polsturm to count roots of a polynomial in an interval using Sturm 
sequences. (This routine requires the polynomial to be squarefree and nonvanishing 
at the left endpoint of the interval; one must write a wrapper function to lift these 
restrictions.) 

• We use the GNU Scientific Library ( GSL ) [9] for computing approximations of roots 
of polynomials over R. 

• We use Cython [B] for migrating some critical routines into compiled C code. 

Additional components we tried out, but did not ultimately use, include Maple (commercial, 
compiled), Numpy (interpreted Python), and Singular (compiled, but its library rootsur 
for Sturm sequences is interpreted at runtime); one side benefit of SAGE is that it makes it 
easy to experiment with many other systems in this manner. Alan Lauder has done some 
further experiments using Magma. 

In this approach, we take the tree in Proposition 13. 21 to be as large as possible, by including 
all tuples not forbidden by (c). This reduces to solving the following problem. 

Problem 4.1. Let R(z) E Q[z] be a polynomial with positive leading coefficient, such that 
R'(z) has all roots real and in [—2^/q, 2^/q). Find all c E Z such that R(z) + c has all roots 
real and in [—2y/q,2y/q\. 

Let x\ < • ■ ■ < Xd-i be the roots of R' counted with multiplicity, and put x = —2^fq 
and Xd = 2^/q. For i — 0, . . . , d, put yi = R(xi). Then the values of c we want are those for 
which 

Z/d- 2l + c>0 (i = 0,...,|d/2j) 
y d -i-2i + c<0 (i = 0,...,|_(d-l)/2j). 

One can interpret this more geometrically by drawing the graph of R over [— 2^/q, 2^fq\. The 
values of c are the negatives of the integral y-values between the highest local minimum and 
the lowest local maximum of R (inclusive), provided that we treat 2^/q as a local maximum, 
and treat —2^fq as a local maximum or minimum depending on whether d is even or odd. 

Our principal method for treating Problem 14. II is to compute numerical approximations to 
the Xi and yi. We throw an exception if these approximations are not sufficiently accurate, 
unless Xi = x i+ i for some i; we can both detect and resolve this case using exact arithmetic. 

Algorithm 4.2. Consider inputs as in Problem l4.ll together with a positive integer p. Using 
GSL, compute numerical approximations x 1 < ■ ■ ■ < x^-i to the roots of R', presumed (but 
not guaranteed) correct to within 2~ p . For i = 1, . . . ,d — 1, put = \_x i 2 p ~ 1 — 1J2~ P+1 
and Si = ri + 2~ p+3 ; also put r = s = —2^/q and r d = s^ = 2^/q. If any of the following 
conditions occur for some i E {1, . . . ,d — 1}: 

• Si > r m ; 

• (-l) d - l R"( ri ) > 0; 

• -R'( r i) an d R'(si) have the same sign; 

then abort or return according as Algorithm 14.41 aborts or returns. If none of the conditions 
occur, put I = — oo and u = +oo. For i — d, d — 1, . . . , in turn: 
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• if n = Si, put t = Rin); 

• if Ti < Si and d — i is even, let i be the value computed by applying Algorithm 14.31 
with [r, s] = [ri, Si] and to = —I, then replace / by max{— t, /}; 

• if Ti < Si and d — i is odd, let £ be the value computed by applying Algorithm 14.31 
with [r, s] = [ri, s^ and t Q = u after replacing R by —R, then replace u by min{t, w}; 

• if now I > u, return the empty set. 

Return the range Z fl [I, «]; this solves Problem 14. II if not aborted. 

Proof. The only thing that needs to be noted here is that failure to invoke Algorithm 14.41 
ensures that the intervals [ri, si), . . . , \r<i-i, Sd-i) are disjoint and contain one root of R' 
apiece, so the input to Algorithm 14.31 is valid. □ 

In order to determine the roundings of the yi, we use exact arithmetic as follows. 

Algorithm 4.3. Let R{z) G Q[z] be a polynomial such that R' has all roots real and distinct. 
Let r, s G Q be such that R"(r) < 0, and the interval [r, s] contains a local maximum of R 
and no other roots of R'. Let t Q G Z U {+oo}. Compute 

t = [R(r)\ 

u = [R(r) + (s-r)R'(r)\. 

If t > to, then return to (this can be checked before computing u). Otherwise, while t ^ u, 
repeat the following: for v = [^p] , if R—v has any roots in [r, s] as determined by polsturm, 
then replace t by v , otherwise replace u by v — 1. Return t; then for a; the unique root of R' 
contained in [r, s], either t > t or t = \_R(x)\ . 

Proof. Since R' has all roots real and distinct, x must be an isolated root of R' . Since x is a 
local maximum for R, R' must undergo a sign crossing at x from positive to negative. Since 
R' has no other roots in [r, s] , R' must be positive in [r, x) and negative in (x, s] . 

The roots of R" interlace those of R' by Rolle's theorem, so in (r, x] we have either zero 
or one root of R" . The root occurs if and only if there is a sign crossing; since R"(x) < 
and R"{r) < 0, we deduce that there is no root, and R"(z) < for all z G (r, x}. 

This implies that R'(r) > R'(z) for z G [r,x]; since R'(r) > 0, 

R{x) = R(r) + / R\z) dz < R(r) + (x - r)R'(r) < R{r) + (s - r)i?'(r). 

This yields the claim. □ 

Note that to a certain extent, taking p small in Algorithm 14. 21 is beneficial to Algorithm I4.3[ 
because it keeps the heights of the rationals r^, Si small. However, it may happen that if p is 
too small, then the gap between the initial values of t and u in Algorithm 14.31 may be quite 
large, and a great deal of time may be wasted narrowing the gap. 

Recall that Algorithm 14.21 does not treat cases of Problem 14.11 in which R! has repeated 
roots, or R\—2yJq)R'{2^fq) = 0; here is a simple treatment. In practice, these cases seems 
to be exceedingly rare; for instance, they do not occur at all in the example of Section [6j 

Algorithm 4.4. Consider inputs as in Problem 14.11 Put T = gcd(i?', (z 2 — Aq)R")\ if T 
is constant, then abort. Otherwise, let S\, . . . , Sk denote the distinct irreducible factors of 
gcd(_R', (z 2 — Aq)R"). Determine whether the quotients upon dividing R by each Si are all 
equal to a single integer — c. If so, use polsturm to check whether R(z) +c has all roots real 
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and in [—2y/q,2y/q\; if so, return the singleton set {c}. In all other cases, return the empty 
set. This solves Problem 14. II if not aborted. 

Proof. Suppose that T is nonconstant and R(z) + c has all roots real and in [—2y/q,2y/q\. 
Let r be a root of T. If r = —2^/q, then by Rolle's theorem, R(z) + c has a root less than 
or equal to —2^/q, hence R(—2^/q) + c = 0. Similarly, if r = 2-Jq, then R(2^/q) + c = 0. If 
—2^/q < r < 2y/q, then r is a root of R" and so must be a multiple root of R'\ by Rolle's 
theorem, r must be a root of R(z) + c. This proves the claim. □ 

5. Second approach: power sums 

Inspection of the enumeration of the maximal tree in some examples suggests that it 
is rather bushy, in the sense of having many vertices with many children but few deep 
descendants. This in turn suggests that a more refined tree construction might be able to 
achieve substantial runtime improvements. Our second approach, implemented in [11] using 
SAGE and components as in the previous section (but again restricted to the case q — 1), 
does this; it is based on estimations of power sums, as in the work of Boyd [I] and subsequent 
authors (most notably [Sj) on searching for polynomials with small Mahler measure. 

Given a polynomial R(z) = ^2™ =0 CiZ l with c n ^ 0, with roots rj, . . . ,r n , the power sums 
of R are defined as 

Sj =r{ + --- + ri (j = 0,l,...). 
They are related to the coefficients of R via the Newton identities: 

j-i 

jc n _j + ^ °n-iSj-i = (j = 1, . . . , n). 

i=0 

In particular, given c n , one can recover c n _i, . . . , c n „j from s\, . . . , Sj via an invertible linear 
transformation over Q. Moreover, the j-th power sum of R(z)+Y^i=o c 'i z% equals Sj—jc' n _j/c n . 
Note that PARI provides a routine polsym to generate the power sums of a polynomial. 

In this tree enumeration, we will generate some nodes which do not actually belong to 
the tree, because they do not satisfy (c); hence our first step when considering a proposed 
node will be to check (c) using polsturm. (Profiling data in some examples suggests that 
this step is a bottleneck in the computation; some improvement may be derived by instead 
using Sturm-Habicht sequences, as described in [2] , or perhaps even using real root isolation 
techniques. We plan to investigate this further.) If (c) is satisfied, and the node is not at 
maximum depth, we enumerate its children by generating and solving an instance of the 
following problem. 

Problem 5.1. Given a polynomial R{z) = Y^7=o ° iZ% w ^ c n 7^ 0, and an integer 1 < j < n, 
find l,u eZ such that for any real numbers d n -j, . . . , do with d n _j 6 Z and R(z) + XT=o ^ iZ% 
having roots in \—2^fq,2^fq\, we have d n _j e [l,u]. 

Note that this problem is somewhat open-ended: if l^u^ is a solution of Problem 15.11 for 
% = 1, . . . , k, then so is I, u for I = maxj{/j},w = minj{wj}. It thus suffices to exhibit a list 
of inequalities satisfied by the coefficients of a polynomial R(z) = XT=o CiZ% with all roots 
in [—2y/q,2^/q\; equivalently, we may exhibit inequalities satisfied by the power sums Si of 
R. Here are some convenient ones; adding additional inequalities should provide even better 
results, although at some point adding a new inequality will eliminate so few cases that it 
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will not be worth the time required to check it. (The linear programming approach in [8] 
may prove helpful in finding good compact systems of inequalities.) 

(1) For % even, 

Si - 4qSi- 2 < (i > 2). 

(2) Let Ti(z) = Y^k=[)ti,k zk be the polynomial of degree i for which T i (2 Ak /gcos6 l ) = 
2 .y/q cos i6 (a rescaled Chebyshev polynomial of the first kind); then 



k=0 



i-2 



k=0 



< 4nqy/q 



(i>0) 



(i > 2). 



(3) Put s[ = ELo (DCV*)* - *** then 



(4) Put < = ELoG)(V?) 



' _fc (— l) fc s fe ; then 



> 



(<>2) 
(i > 1). 

(z>2) 
fi > 1). 



6. An example 

Here is an example of the basic problem, excerpted from [H §4.2], and some results obtained 
using the algorithms and implementations described above. 

Consider the smooth quartic surface X in the projective space over the finite field F3 
defined by the homogeneous polynomial 

x 4 — xy 3 + xy 2 w + xyzw + xyw 2 — xzw 2 + y 4 + y 3 w — y 2 zw + z 4 + w 4 . 

(As described in [TJ, this polynomial was chosen essentially at random except for a skew 
towards sparseness.) Since X is a K3 surface, the Hodge diamond of X is 

1 


1 20 1 

1 

and the Hodge polygon of primitive middle cohomology has vertices (0, 0), (1, 0), (20, 19), (21, 21). 
Consequently, the zeta function of X has the form 



rpf. 



Cx(T)=exp V— #X(F 3 n 
\ ^— 1 n 



1 



v n=l 



1 -T)(l -3T)(1 -9T)R(T) 



where R(T) G Z[T] is a polynomial of degree 21 such that R(0) = 1, the complex roots of 
R lie on the circle |T| = 3 _1 , and (by an inequality of Mazur) the Newton polygon of R lies 
above the Hodge polygon. In particular, the polynomial S(T) = 3R(T/3) is root-unitary 
and has integral coefficients. 
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Define 

So CD = 3T 21 + 5T 20 + 6T 19 + 7T 18 + 5T 17 + 4T 16 + 2T 15 - T 14 - 3T 13 - 5T 12 

- 5T 11 - 5T 10 - 5T 9 - 3T 8 — T 7 + 2T 6 + 4T 5 + 5T 4 + 7T 3 + 6T 2 + 5T + 3; 

one easily checks that So is root-unitary. By explicitly enumerating X(¥ q n) for n < 5, one 
finds 

S(T) = S (T) (rnodT 6 ); 
by a 3-adic cohomology computation described in pQ, one finds 

5(T) = Sq(T) (mod 3 5 ). 

Having performed these computations, one wants to verify whether these restrictions suffice 
to ensure S(T) = Sq(T). Moreover, one also wants to know to what extent they can be 
weakened while still forcing S(T) = Sq(T), as the enumeration and cohomology calculations 
become significantly more cumbersome as the strength of their results is forced to increase. 

The result obtained here is that already the conditions that S(T) G Z[T], S(T) is root- 
unitary, and 

S(T) = S Q (T) (mod 3 2 T X ) 
force S = So. Note that already the congruence S(T) = Sq(T) (mod 3) implies that S must 
be reciprocal rather than antireciprocal, so we may as well put 

P(T) = S(T)/(T + 1), P (T) = S (T)/(T + 1). 

Then the conditions we are interested in are that P{T) e Z[T], P(T) is root-unitary and 
reciprocal, and 

(6.1) P{T) = P (T) (mod ??T j ) 

for various The asserted result is that these conditions for i — 2,j — 1 force P = Pq. 
(It turns out that % = 1 does not suffice even with j = 10.) 

We checked the sufficiency of the conditions for i = 2,3,4,5 and j = 1,2,3,4,5 by run- 
ning the implementation [T2] on one Opteron 246 CPU (64-bit, 2 GHz) of the computer 
dwork.mit . edu. The machine has 2GB of RAM. 

The timings and sizes of the computations for various initial constraints are summarized 
in Table [1] (using root-finding, as in Section HJ with rounding precision p = 32) and Table [2] 
(using power sums, as in Section [5]). Each entry consists of the number of CPU seconds 
for the calculation, rounded up to the nearest tenth of a second, followed by the number of 
leaves (terminal nodes) in the tree over which we exhausted. Note the significant savings 
achieved by the second approach. (We did some additional experiments combining the two 
approaches, but the power sum method by itself seemed to outperform hybrid methods.) 

TABLE 1. Timings for recovery of P given (I6.ip . using root-finding. 





3 2 


3 3 


3 4 


3 b 


T i 


564.2/1011788 


2.2/3858 


0.1/38 


0.1/2 


T 2 


267.9/501620 


2.2/3784 


0.1/38 


0.1/2 


T 3 


4.3/4714 


0.1/63 


0.1/6 


0.1/1 




1.8/1838 


0.1/51 


0.1/6 


0.1/1 


y5 


0.7/612 


0.1/32 


0.1/5 


0.1/1 
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Table 2. Timings for recovery of P given (I6.ip . using power sums. 





3 2 


3 3 


3 4 


3 b 


T 1 


1.9/1157 


0.1/6 


0.1/1 


0.1/1 


T 2 


0.6/347 


0.1/2 


0.1/1 


0.1/1 


T 3 


0.3/117 


0.1/1 


0.1/1 


0.1/1 




0.2/53 


0.1/1 


0.1/1 


0.1/1 


j~>5 


0.1/23 


0.1/1 


0.1/1 


0.1/1 



7. Another example 

The previous example shows the superiority of the power sum method over the root-finding 
method. This suggests trying a larger example to test the limits of the power sum method; 
here is an example provided by Alan Lauder. 

The reciprocal polynomial 

P(T) = 2401T 56 - 343T 55 - 54 39T 54 - 1050T 53 + 7156T 52 + 50 43T 51 - 5829T 50 - 7990T 49 

+1437T 48 + 6348T 47 + 2115T 46 - 332T 45 - 1756T 44 - 4639T 43 - 1802T 42 + 39 38T 41 

+4762T 40 + 16T 39 - 3366T 38 - 2658T 37 - 2051T 36 + 1572T 35 + 58 1 0T 34 + 2097T 33 

-5558T 32 - 3955T 31 + 2598T 30 + 1931T 29 - 831T 28 + 1931T 27 + • • • 

is root-unitary; it arises from a 7-adic cohomology calculation of the primitive middle coho- 
mology of an elliptic surface over F7 with Hodge diamond 

1 


4 49 4. 

1 

As in the previous example, we ask whether a reciprocal root-unitary polynomial Pq{T) 
satisfying 

(7.1) P(T) = P (T) (mod 7*T J ) 

necessarily equals P(T). In the following, each expression (A/B) indicates that the indicated 
computation required A CPU seconds and encountered B terminal nodes. 

• For i = 2, P = P is not forced for j = 28 (0.7/15). 

• For z = 3, P = P is forced for j = 25 (336.2/355435) but not for j = 24 
(711.7/755544). 

• For i = 4, P = P is forced for j = 16 (331304.0/196405710). We were unable to 
find any value of j for which Pq = P is not forced. (For comparison, the complexi- 
ties for j = 17, 18, 19,20 were 61787.7/36665858, 12464.5/7334642, 2275.5/1349860, 
392.6/232783.) 

• For i = 5, P = P is forced for j = 1 (93.7/13513). 

One can explain this behavior heuristically by imposing only the condition that (in the 

notation of Section[5]) | X]l=o^ fc,Sfc l — ^ n ^ or * — 0- This restriction constrains the coefficient 
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of T- 5 to a range of size (4 • 28)/(7 4 • j). Once i is big enough that this range typically takes 
only a few elements, we can expect to be able to force Pq = P. 

8. The case of nonsquare q 

As noted earlier, our implementations so far have only covered the case q = 1. It is easy 
to reduce to this case from any case in which q is a square. For q not a square, there are 
several ways to proceed; we do not know which of these is best. 

• One may repeat the methods as written above, but using exact arithmetic in the 
quadratic field Q(y/q). 

• One may replace yfq by an upper approximation by a rational number s and look for 
polynomials with roots in [—2s, 2s], then screen out those which do not have roots 
in [-2^,2^]. 

• One may consider the polynomial S defined by S(z 2 ) = R(z)R(—z). 

9. Further comments 

Note that SAGE runs primarily in the interpreted language Python, although many of 
its components either are compiled libraries, or have been migrated to C using Cython. 
(Indeed, the latter progress has been ongoing, and this can be detected in the runtimes 
of our algorithms under different versions of SAGE.) It is thus fair to ask whether some 
additional optimization could be achieved by porting everything to a compiled environment. 
We have already built in some savings by performing most polynomial manipulation in PARI 
with limited conversions to/from SAGE, and by porting some key subroutines into Cython; 
it is not clear how much more room there is for improvement on this front. 

Our depth-first search is implemented using a queue rather than recursion; this has the 
advantage of making it easily amenable to parallelization. Although quite sophisticated 
strategies have been devised for scheduling in the context of tree traversal (e.g., [17]), even 
implementing some simple scheduling mechanisms, such as work-stealing, would be helpful 
in a multiprocessor environment. Starting with version 2.0, SAGE includes a subsystem 
called DSAGE (Distributed SAGE), which may facilitate this sort of simple parallelization. 

One can use similar search techniques for polynomials with roots bounded in a convex 
subset of the complex plane, since the Gauss-Lucas theorem asserts that this property is 
also preserved by taking derivatives. We have not experimented with this in any detail. 
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